function plotrocks(filenumber)
close all;
cmap=[
    1.0000    1.0000    1.0000
    0.9167    0.9167    1.0000
    0.8333    0.8333    1.0000
    0.7500    0.7500    1.0000
    0.6667    0.6667    1.0000
    0.5833    0.5833    1.0000
    0.5000    0.5000    1.0000
    0.4167    0.4167    1.0000
    0.3333    0.3333    1.0000
    0.2500    0.2500    1.0000
    0.1667    0.1667    1.0000
    0.0833    0.0833    1.0000
         0         0    1.0000
         0    0.0667    1.0000
         0    0.1333    1.0000
         0    0.2000    1.0000
         0    0.2667    1.0000
         0    0.3333    1.0000
         0    0.4000    1.0000
         0    0.4667    1.0000
         0    0.5333    1.0000
         0    0.6000    1.0000
         0    0.6667    1.0000
         0    0.7333    1.0000
         0    0.8000    1.0000
         0    0.8667    1.0000
         0    0.9333    1.0000
         0    1.0000    1.0000
    0.0588    1.0000    0.9412
    0.1176    1.0000    0.8824
    0.1765    1.0000    0.8235
    0.2353    1.0000    0.7647
    0.2941    1.0000    0.7059
    0.3529    1.0000    0.6471
    0.4118    1.0000    0.5882
    0.4706    1.0000    0.5294
    0.5294    1.0000    0.4706
    0.5882    1.0000    0.4118
    0.6471    1.0000    0.3529
    0.7059    1.0000    0.2941
    0.7647    1.0000    0.2353
    0.8235    1.0000    0.1765
    0.8824    1.0000    0.1176
    0.9412    1.0000    0.0588
    1.0000    1.0000         0
    1.0000    0.9444         0
    1.0000    0.8889         0
    1.0000    0.8333         0
    1.0000    0.7778         0
    1.0000    0.7222         0
    1.0000    0.6667         0
    1.0000    0.6111         0
    1.0000    0.5556         0
    1.0000    0.5000         0
    1.0000    0.4444         0
    1.0000    0.3889         0
    1.0000    0.3333         0
    1.0000    0.2778         0
    1.0000    0.2222         0
    1.0000    0.1667         0
    1.0000    0.1111         0
    1.0000    0.0556         0
    1.0000         0         0
    0.5000         0         0 ];

fullscreen = get(0,'ScreenSize');
hf=figure('Position',[0 -50 fullscreen(3) fullscreen(4)]);
xzoomlim=1e6*[0.8930    1.3046];
yzoomlim=1e3*[-0.0523    1.0468];

width=10;
%center=1060; xwidth=500e3; ylim=[0 200e3];
% papersizepos=[0 0 width width*(ylim(2-ylim(1)))/xwidth];
% set(hf,'unit','centimeter','position',papersizepos,'paperunits','centimeter','paperposition',papersizepos);
% h1=axes('fontsize',4,'units','centimeter','pos',papersizepos);


load('..\..\crosstopodata2.mat');
N=readdata(0,1);
%     S=readdata(0,2);
%     T=readdata(0,3);
%     multigridlevels=readdata(0,4);
%     bound=readdata(0,5);
time=readdata(0,6);
%     dtnom=readdata(0,7);
dx=readdata(0,8);
dz=readdata(0,9);
mx=readdata(0,10);
mz=readdata(0,11);
%     mcs=readdata(0,12);
%     mct=readdata(0,13);
%     mtype=readdata(0,14);
%     mTemp=readdata(0,15);
%     mstressxx=readdata(0,16);
%     mstressxz=readdata(0,17);
%     mstrainxx=readdata(0,18);
%     mstrainxz=readdata(0,19);
%     mfinitestrain=readdata(0,20);
%     vx=readdata(0,21);
%     vz=readdata(0,22);
%     P=readdata(0,23);
Nsurfaces=readdata(0,24);
surfaces=readdata(0,25);

mcs=zeros(1,N); mct=zeros(1,N);
xdivision=1; xplotcellsize=sum(dx)/xdivision;
x=0; for s=0:xdivision-1 ind=find(mx>=x & mx<x+xplotcellsize); mcs(ind)=s; x=x+xplotcellsize; end
zdivision=150; zplotcellsize=sum(dz)/zdivision;
z=0; for t=0:zdivision-1 ind=find(mz>=z & mz<z+zplotcellsize); mct(ind)=t; z=z+zplotcellsize; end
%Ia=int32(find((-1).^mcs.*(-1).^mct>0));
Ib=int32(find((-1).^mcs.*(-1).^mct<=0));
clear('mcs','mct');

%mtotalstrain=zeros(N,1); 
oldtime=time;

surfaces=readdata(1,25);
reftopo=surfaces(1,3)-surfaces(:,4);

for count=filenumber

    [xplates,gz]=calcgrav(count,1000,0);
    
    N=readdata(count,1);
    S=readdata(count,2);
    T=readdata(count,3);
    %     multigridlevels=readdata(count,4);
    %     bound=readdata(count,5);
    time=readdata(count,6);
    %     dtnom=readdata(count,7);
    dx=readdata(count,8);
    dz=readdata(count,9);
    mx=readdata(count,10);
    mz=readdata(count,11);
    %     mcs=readdata(count,12);
    %     mct=readdata(count,13);
    mtype=readdata(count,14);
    %mTemp=readdata(count,15);
    %mstressxx=readdata(count,16);
    %mstressxz=readdata(count,17);
    %mstrainxx=readdata(count,18);
    %mstrainxz=readdata(count,19);
    %mfinitestrain=readdata(count,20);
    %     vx=readdata(count,21);
    %     vz=readdata(count,22);
    %P=readdata(count,23);
    Nsurfaces=readdata(count,24);
    surfaces=readdata(count,25);
    
    %center=1060e3;
    dt=time-oldtime; oldtime=time;
    x=zeros(1,S); z=zeros(1,T);
    for s=1:S-1 x(s+1)=x(s)+dx(s); end
    for t=1:T-1 z(t+1)=z(t)+dz(t); end
    figure(1); clf;
    
    sealevel=surfaces(1,3);
    ind=find(x<0.6e6); surfaces(ind,Nsurfaces)=surfaces(ind,4); %rettelse
    clear ind;

    surfacehandle=subplot(4,1,1);
    topography=sealevel-surfaces(:,Nsurfaces);

    [ax,h1,h2]=plotyy(x/1e3,topography+773,xplates/1e3,gz/1e-5); hold(ax(1)); hold(ax(2)); set(h1,'color','g','linewidth',2); set(h2,'color','b','linewidth',2);

    %     plot(x/1e3,topography+773,'g','linewidth',1); hold on;
    plot(ax(1),profilex+630,profile_t,'g--','linewidth',1); plot(ax(2),profilex+630,profile_g,'b--','linewidth',1);
    set(ax(1),'ycolor','k'); ylabel(ax(1),'Topography (m)'); set(ax(1),'xtick',[]); set(ax(1),'ylim',[0 7000]); set(ax(1),'ytick',0:1000:7000); set(ax(1),'box','off');
    set(ax(2),'ycolor','k'); ylabel(ax(2),'Free air gravity (mgal)'); set(ax(2),'xtick',[]); set(ax(2),'ylim',[-400 400]); set(ax(2),'ytick',-400:200:400); set(ax(2),'box','off');
    %
    
    mTemp=readdata(count,15);
    [grid,xgrid,zgrid]=bilingrid(mTemp,mx,mz,dx,dz,1,1);
    contvals=1200:20:1800; cont=contourc(xgrid,zgrid,grid',contvals);
    clear('mTemp');
    avtemp=0*zgrid;
    for t=1:length(avtemp);
        avtemp(t)=mean(grid(:,t));
    end
    ztemp=zgrid;
    clear('grid');
    
    mtype(find(mtype==6))=1; mtype(find(mtype==0))=-1; %mtype(find(mtype==3))=2; mtype(find(mtype==4))=3; mtype(find(mtype==5))=4;
    mtype(Ib)=mtype(Ib)+0.5; mtype(mtype<0)=-1; %mtype(mtype>5.1)=5; 
    [grid,xgrid,zgrid]=bilingrid(mtype,mx,mz,dx,dz,2,2); %grid(find(grid<0.75))=nan;
    clear('mx','mz','mtype','Ib');
    grid=round(2*grid);
    [Zgrid,Xgrid] = meshgrid(zgrid,xgrid);
%     xrange=find(xgrid>center-250e3 & xgrid<center+250e3); zrange=find(zgrid>0 & zgrid<200e3);
%     Xgrid=Xgrid(xrange,zrange); Zgrid=Zgrid(xrange,zrange); grid=grid(xrange,zrange);

    %grid(1:5,1:5)=5.5;
    lithohandle=subplot(4,1,2:4);
    h(1)=pcolorim(gca,Xgrid'/1e3,(Zgrid'-sealevel)/1e3,grid',1,1); set(gca,'ydir','rev'); hold on; clear('Xgrid','Zgrid','grid');
    xlabel('Distance (km)'); ylabel('Depth (km)');
    colormap(cmap); %caxis([-1 11]);
    %hc=colorbar; oldcbarpos=get(hc,'pos');
    hold on;

    pos=1; allypoints=cell(1,length(contvals));
    while pos<length(cont)
        step=cont(2,pos);
        i=find(contvals==cont(1,pos));
        xcont=cont(1,pos+1:pos+step); ycont=cont(2,pos+1:pos+step);
        %plot(xcont,ycont-sealevel,'color',[0.99 0.99 0.99],'linewidth',0.1);
        plot(xcont/1e3,(ycont-sealevel)/1e3,'k','linewidth',0.1);
        pos=pos+step+1;
        allypoints(i)={[cell2mat(allypoints(i)) ycont]};
    end
    xlim=[0 max(x)/1e3]; ylim=[-20 641]; ratio=(ylim(2)-ylim(1))/(xlim(2)-xlim(1));
    axis([xlim ylim]);
    set(gca,'box','off');
    set(ax(1),'xlim',get(gca,'xlim')); set(ax(2),'xlim',get(gca,'xlim'));
    lithopos=get(lithohandle,'pos'); lithopos(4)=lithopos(3)*ratio; set(lithohandle,'pos',lithopos);
    
    newxlim=[900 1200]; newylim=[-10 250];
    fullscreen = get(0,'ScreenSize'); figure('Position',[0 -50 fullscreen(3) fullscreen(4)])
    cdata=get(h(1),'cdata'); imagesc([x(1) x(end)]/1e3,([z(1) z(end)]-sealevel)/1e3,cdata); axis equal; axis([newxlim newylim]); colormap(cmap); %caxis([-1 11]); set(gca,'box','off');
    print -depsc -painters fig1b; close(gcf);
    
    figure(hf);
    newaxespos=get(lithohandle,'pos');
    newaxespos(3)=newaxespos(3)/4;
    axes('pos',newaxespos);
    plot(avtemp,(ztemp-sealevel)/1e3,'r','linewidth',2);
    axis([0 1900 ylim]);  set(gca, 'color', 'none'); set(gca,'ydir','rev'); set(gca,'ytick',[]); set(gca,'XAxisLocation','top');  set(gca,'box','off'); xlabel('Horisontally averaged temperature (^{\circ}C)');
    

    
    xp=[newxlim(1) newxlim(1) newxlim(2) newxlim(2) newxlim(1)]; zp=[newylim(1) newylim(2) newylim(2) newylim(1) newylim(1)];
    %zp=zp+500;
    plot(lithohandle,xp,zp,'k','linewidth',0.3);
    
    newaxespos=get(lithohandle,'pos'); rightmargin=newaxespos(1)+newaxespos(3);
    newaxespos(3)=newaxespos(3)/4; newaxespos(1)=rightmargin-newaxespos(3);
    mx=readdata(count,10); mz=readdata(count,11);
    mstressxx=readdata(count,16); mstressxz=readdata(count,17); mstresstot=(mstressxx.^2+mstressxz.^2).^0.5; clear('mstressxx','mstressxz');
    [grid,xgrid,zgrid]=bilingrid(mstresstot,mx,mz,dx,dz,2,2);
    avstress=0*zgrid;
    for t=1:length(avstress);
        avstress(t)=mean(grid(:,t));
    end
    [nothing stresspos]=min(abs(xgrid-1.25e6));
    oldaxespos=get(lithohandle,'pos');
    newaxespos(1)=oldaxespos(1)+xgrid(stresspos)/xlim(2)/1e3*oldaxespos(3);
    axes('pos',newaxespos); plot(grid(stresspos,:)/1e6,(zgrid-sealevel)/1e3,'r','linewidth',2); axis([0 1000 ylim]);
    set(gca, 'color', 'none'); set(gca,'ydir','rev');  set(gca,'XAxisLocation','top');  set(gca,'box','off'); set(gca,'ytick',[]); xlabel('\sigma_{II}-profile (MPa)'); %set(gca,'YAxisLocation','right'); set(gca,'xdir','rev');
    
    
    set(gcf,'position',[0          -7        1115        1130]);
    set(gcf,'paperPositionMode','auto');
    %set(gcf,'paperposition',[0          -7        1115        1130]);
    print -depsc -painters fig1a;
    


end
end